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ABSTRACT 

We present a unified parameterization of the fitting functions suitable for density profiles 
of dark matter haloes or elliptical galaxies. A notable feature is that the classical Einasto 
profile appears naturally as the continuous limiting case of the cored subfamily amongst the 
double power-law profiles of IZhaoL Based on this, we also argue that there is basically no 
qualit ative difference between halo models well-fitted by the Einasto profile and the standard 
iNFWl model. This may even be the case quantitatively unless the resolutions of simulations 
and the precisions of fittings are sufficiently high to make meaningful distinction possible. 

Key words: dark matters - method: analytical. 



1 INTRODUCTION 

Although it is often claimed that the density profiles of dark mat- 
ter haloes found in simulations are well described by simple dou- 
ble po wer-law functions of radii (see e.g.. lNavarro. Frenk & White] 
1 19951 henceforth NFW), these are usually based on the fitting with 
limited radial coverage. A possible manifestation of this limitation 
is that these functions, when they are extrapolated to the centre, typ- 
ically possess a singularity, which is generally unphysical. Real dy- 
namical systems almost certainly have finite phase-space densities, 
finite number densities, and finite escape speeds even at the cen- 
tre, as argued by fundamental physics: (1) any halo with fermions 
are limited in phase-space density by the Pauli exclusion principle; 
(2) any annihilating dark matter particle is limited in the number 
density by the annihilation cross section; (3) any stellar system is 
limited by its relaxation time; and (4) the escape speed is limited at 
least by the speed of light. To ameliorate the singularity whilst in- 
troducing the fewest extra parameters is to make the density profile 
more flexible in its parameterization. 

In fact, the results from mo re recent high resolution simulations 
(see e.g.. iNavarro et al] l2010h appear to suggest that a di fferent 
class of fitting formula such as that o f lEinastol JT965LI 1969h intro- 
duced earlier for spherical stellar systems may be more close to the 
'reality'. The Einasto profiles by construction have finite densities, 
but they are less user friendly given their exponential behaviours. 
This motivates us to look for smooth transitions from double power 
laws to the Einasto models. Although they seem to be qualitatively 
different at first glance, the fitting formulae for the double power 
laws and the Einasto profiles are indeed closely related. For exam- 
ple, the Gaussian and the Plummer profile are limiting cases of the 
so-called beta profile for hot X-ray gas haloes, oc (1 + r 2 / a 2 \-(P+ 1 / 2 ) 
for f3 — > oo or /3 = 2 respectively. This hints to us that it may also 
be possible to connect the entire family of double power laws and 



the Einasto profiles. We shall show that this is achieved by simply 
redefining parameters of the double power-law models. This is also 
analogous to relating an isothermal sphere to polytropic models as 
in the infinite-index power law reducing to an exponential. 

We argue that this connection is more than a mere mathemati- 
cal trick, for it suggests that a cored double power-law profile with 
a large outer power-index value and the Einasto profile should be 
qualitatively indistinguishable if the fitting is mostly weighted near 
and around the central regions. Moreover, the resolution limit in- 
dicates that the fittings to dark halo profiles, be they numerical 
simulations or some observational proxies, are mostly concerned 
with the behaviour around the 'scale radius', that is, correspond- 
ing to the transition regions in the model. Hence, combined with 
the uncertainties in the fittings, it is expected that most cusped dou- 
ble po wer-law profile s with a modera te outer slope, such as the 
NFW profile and the iHernauistl jl99fj|) model, may not be well- 
distinguished from some cored profiles with an extreme value of 
the outer slope or even exponential fall-off if the variations of the 
logarithmic slope in the latter models are made to be sufficiently 
slower than the former. Here we shall show that, in the context of 
our extended family of models, this is indeed the case. 

Despite widespread use of the double power-law profile, this pos- 
sibility has not been investigated in detail, for in practice most fo- 
cus on the double power law has been limited to the two limiting 
power indices and scarce attention has been paid to the parame- 
ter controlling the sharpness of the transitio n. The connection of 
the double power-law profiles to the lEinasto profile, for which the 
variation of the shape parameter has received more attention, thus 
also highlights the role of the corresponding parameter in the dou- 
ble power law, and subsequently the above-mentioned difficulty of 
distinguishing fitting functions. 
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Figure 1. Density profiles of selected members of the family. All models 
drawn have a = 1 and p = 1 but the values of the parameter s vary with 
c5 = 3,4,5,9, 15 and s = 0, where c5 = s . The model s with (a, 6, p) = 
(1,3, 1) a nd (a, S, p) = (1, 4, 1) correspond to the so-called NFW profile and 
iHernquisl model, respectively. In the top panel, the plots are normalized for 
a fixed g s and r s whilst they are for Q-2 and r_2 in the bottom panel. In 
addition, in the bottom panel, they are further scaled by r~ 2 so that the 
horizontal tangent in the figure indicates that local power index of -2. 



2 FITTING FORMULAE FOR DENSITY PROFILES 
2.1 The Double power law 

A widely-used fitting formula for the density profiles of the dark 
haloes found in cosmological sim ulations is of the form of a ' bro- 
ken' or double power-law profile (Hernquist 199(il: IZhadll996h : 



diogp aa p + 6r p 



dlog r 



a p + rP 



(1) 



Here, the parameters a < 3 and c5 correspond to the power indices 
of the density profile at the centre and towards the asymptotic limit. 
By considering only such profiles with vanishing density at infinity, 
with possibility of escape (i.e., the potential at infinity can be set to 
a finite value), or with a finite total mass, allowed values of the 
parameter c5 may be restricted to be c5 > 0, 6 > 2, or S > 3. The 
parameter p > controls the width of the transition region between 
the two limiting behaviours; that is, the transition becomes sharper 
as p gets larger. Finally, the scale length a in equation (TJ specifies 
the radius at which the local logarithmic density slope is the same 
as the arithmetic mean of the two limiting values. 

The well-known specia l case s of these models include: the 
ISchusted dl884h4Plummeildl91lh s phere (a = 0, 6 = 5, p = 2: an 
ind ex-5 polytrope); t he Jaffe (1983) model (a = 2, S = 4, p = 1); 
the iHeniquistl dl99d) model (a = 1, S = 4, p = 1); the lNFWl 
pr ofile (a = 1, 6 = 3 , p = 1); and the isotropic analytic solution 
of lAustin et all J2005t) and lDehnen & McLaughlin! d2005l) (or = |, 
6 = y, p = a)- Also notable are t he families of models : the y- 
sphere studied by Dehnen ( 1 9931) and lTremaine et alJdl994h (c5 = 4, 
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Figure 2. Same as Fig. \T\ but with different members of the family with 
varying p. Again all models shown have r~' cusp (a = 1) but only those 
with two different outer slopes, 6 = 3 (drawn in thick lines) and an expo- 
nential fall-off (j = 0; thin lines) are chosen. Varying p values are indicated 



by different line types: dashed (p 
and dotted (p = | ). 



2), solid (p = 1), dot-dashed (p = ±), 



p = 1); the generalized NFW profil es ofiNavarro etaD ( ll996l . ll997h 
(6 = 3, p = 1); the ^-sphere of lzhad dl996h (a = 1, p = 1: 
which corresponds to alte rnative generalized NFW profiles stud- 
ied by lEvans & Z3I2006); the hypervirial family of lEvans"& a3 
(2005) (a = 2 - p, 6 = p + 3: which was originally i ntroduced 
as the generalized isochronous mod el bvlVeltmannlll979l): and the 
phase-space power-law solutions of Dehnen & McLaughlin ( 2005) 
(a = 1 - | + ySo, S = 3 + p: where /? is the anisotropy parameter at 
the centre). 

2.2 The Einasto profile 

iNavarro et al. ||200i) introduced a different class of fitting formulae 
for the density profile of dark haloes, namely, 

dlog£> 



pocexpi-!? — 



dlogr 



-pn— ■ 

a p 



(2) 



Here p is the parameter controlling how rapidly the logarithmic 
density slope varies. Note the p = 2 case corresponds to the 
isotropic Gaussian. The choices of the scale length a and the con- 
stant t] are not independent but the physical definition of one spec- 
ifies the other. For example, if a = r_ 2 where r_ 2 is the radius at 
which the logarithmic density slope is '-2', then r\ = Th is fam- 
ily^ usually referred to as the Einasto profile after Einastc] \ 19651 
Il969h. It also ha s th e same form as t he fitt ing function proposed 
bv lSersid d 19681) and lde Vaucouleurj d 19481) (the latter is a partic- 
ular case of the former), which have been used to fit the surface 
bri ghtness profile of ell iptical (or spheroidal components of) galax- 
ies dMerritt et al. 2 0051) . 
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Figure 3. Rotation cu rve for sel ect models; InfwI profile (a = 1, 5 = 3, 
p = 1 ; solid lines), the lHernquisJ model (a = 2, S = 4, p = 1; dashed lines), 
Schuster- Plummet sphere (a = 0, 6 = 5, p = 2; dotted lines), r -1 cusped 
exponen tial fall-off (a = 1, * = 0, p = 1; dot-dashed lines), and Sersic- 
lEinastd profile with p = I (a = 0, s = 0, p = j). The plots are normalized 
to have th e same circular speed at r = r_2 ■ Among the models shown, only 
the InfwI profile has an infinite total mass. All remaining models have a 
finite total mass and so their circular speed speeds fall off as 1 / Vr. 



3 AN EXTENSION OF THE DOUBLE POWER LAW 

Functions used to fit the logarithmic density slope in equations Q3 
and are first-order rational functions of r p . Here we consider a 
slight variant version of equation given by 

dlogg _ arj! + r p 
d log r rf + srP 

This is actually the most general first-order monomial rational func- 
tion of r p encompassing both equations and that may be used 
to fit any physical profile of the logarithm density slope. Here we 
assume r s > without any loss of generality (the r s = case cor- 
responds to a pure power law, which may be described by other 
parameter combinations). This reduces to equation through the 
parameter transformation 6 = s" 1 and a 1 ' = Sr P , but it now in- 
cludes the models with s = (i.e., S = oo) in places of those 
with (5 = 0, which would represent unphysical asymptotically non- 
vanishing density profiles. The (a, s) = (0, 0) models correspond to 
those of equation with a p = prjr p . That is to say, equation 
with a = (a cored double power law) reduces to equation by 
taking the limit of S — > oo and a — > oo whilst maintaining a p /6 = rf 
constant (and redefining a p = prjr p afterwards). 

In equation 0, the parameter a < 3 still specifies the central 
power index whereas s > is now the reciprocal of the asymp- 
totic index. The limits for the finite total mass, the escapability, and 
the well-defined integrated density along any line of sight are given 
respectively by < s < | , < s < | , and < s < 1 . We may intro- 
duce a further restriction s < a~ [ by limiting ourselves to models 
whose outer logarithmic slope is steeper than the inner value. The 
parameter p > again controls the steepness of the index variation 
and the breadth of the transition region. On the other hand, although 
the parameter r s still specifies the nominal scale length, it lacks any 



immediate physical interpretation. By allowing an infinite asymp- 
totic power index, the arithmetic mean of the limiting indices is not 
necessarily defined (note the scale length defined as eq. Q~]becomes 
infinite). An alternative physical scale length may be defined such 
as the radius at which the local logarithmic slope is the same as a 
given fixed value, e.g., 



r - o r, . 



(4) 

1 — ns 6 — n 
Here the corresponding logarithmic slope is {-n) < 0. This is de- 
fined for the models for a < n and s < rT l , and simplifies, if s = 0, 
to r? M = (n - a)r{. The most common choice for the reference 
logarithmic slope is (-n) = (-2), which corresponds to that of a 
singular isothermal sphere. 



3.1 Density profiles 

The actual functional form of the density profile is obtained by in- 
tegrating equation 0. If s > 0, this results in 



6-a K 

e(r) = £> s (1 + s) o — 



1 + 



Sri 



(5a) 



1 + 



2-a 
6-2 



1-2 



which reproduces the double power-law profile (we have used the 
substitution 5 = s~ l > for clarity). The integration also introduces 
an additional constant of integration, namely, the scale constant. We 
have chosen this scale constant in reference to the scale length, that 
is, e(r s ) = q s and p(r_ 2 ) = Q-2- The relation between them is 

^ = o^/i+*r (5b ) 

q s ) (2-a) a \S-aj 

We also use q s = (1 + sY 6 ~ a ^ /p Q s , which is useful for compactly writ- 
ing down many properties of haloes. This is also related to the pro- 
portionality constant of the usual expression of the double power 
law in equation such that C = ga" = g s r" given a p = dr{ . 
With s = 0, the integration leads to the density profile given by 



Q(r) 



:£> s e 1/p -J-exp 



Q-2 



exp 



2 



P 



— I -1 

r-2 



and the scale constants related by 
gj (2 - ay ' 



(6a) 



(6b) 



We also define j5 s = e ' p g s , which is the limit of the prior definition 
as s — > 0. The transition of models from s > to s = is also 
smooth - equation J6ab is indeed the limit of equation < f5at as s — > 
(S — > oo). The profile in equation < f6ab is a generalization of the 
Einasto profile of equation allowing a central density cusp. The 
classical profile consists of the cored members (a = 0) of the family 
with the index n = p~ l whilst a cusped model (a = -1.85) has b een 
recently used to model the Galactic spheroid JWang et alj|2012l) . 

Figures [T] and [2] present the log-log plots of density profiles for 
select members of the family. Here the shape of the profiles is fixed 
by three of the parameters (a, s = 6~ [ ,p), whereas varying the scale 
length or constant simply translates the graph horizontally (for the 
scale length) or vertically (for the scale constant). Figure[T]demon- 
strates that the transition of an extreme power law to an exponential 
fall-off is indeed natural. The role of parameter p as characterizing 
the breath of the transition region is highlighted in Figure [2] 
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Figure 4. Energy part of the distribution function (upper panel) and the differential energy distribution (lower panel) of the constant-anisotropy Hernauist model 
with a black hole for (3 = j (left) and /3 = - i (right): thin solid lines (no black hole), thick dotted lines (M, = M tot /100), thick solid lines (M. = M tot /10), 
thin dotted lines (M. = M tot /2), short-dashed lines (M. = M tot ), long-dashed lines (M. = 3M tot /2), dot-dashed lines (M. = 2M tot ) where M tot is the total 
dark halo mass. All models are normalized to the same M tot . Note < £ < GM tot /a when there is no black hole whereas £ > otherwise. 



3.2 Basic properties 

Deriving basic properties of these models such as the enclosed mass 
profiles M r , the rotation curves (i.e., the circular speed) v c , the po- 
tential if/ etc. is basically an exercise in integration. These integrals 
are trivial to evaluate numerically and usually result in special func- 
tions of radii no more complicated than beta and gamma functions. 
We refer Appendix [A] for more systematic studies, which may be 
of some peripheral interest. Here instead we simply provide the ro- 
tation curves for some familiar models and those with exponential 
fall-off in Figure [3] Since the total mass is finite for 6 > 3, the 
circular speeds for all plotted models fall off as Keplerian - i.e., 
v 2 e f?(r~'|3 - for large enough radii, except for thelNFW profile, 
whose mass profile grows logarithmically and thus v 2 falls off like 
r~ l In r. For models with S < 3, it is obvious that the circular speed 
would behave like v 2 e 0(r 2 ~ 6 ) as r — > <x>. None of the plotted mod- 
els possesses a cusp diverging faster or as fast as that of a singular 
isothermal sphere. Hence v 2 . e 0{r 2 ~ a ) and so v 2 — » as r — > in 
Figure [3] If a < 2 < S like those plotted, v c attains its maximum 
at the radius corresponding to the solution of M r = Anr'gir) (i.e., 
d log M,-/d log r = 1). This occurs around r « r_ 2 for the models 
plotted in Figure [3] Either if 5 is not so much different from the 
nominal value of '-2' or if the profile has an extended transition 
region (achieved by small p), the variation of v c then become suffi- 
ciently slow so that the behaviour of the circular speed can mimic 
the flat rotation curve near the region around r » r_2- 

We also find that M r 6 <9(r 3 ~") as r — * for a < 3, which is 
obvious from the fact that dM, /dr = 4nr 2 g e &(r 2 ~") as r — > 
(assuming the null integration constant). Realistically, any model 
with 2 < a < 3 cannot be strictly physical extending down all the 
way to the centre, for the corresponding M r e 0{r 3 ~ a ) would im- 



1 f(x) s @[g(x)] as x — > xq if and only if lim r _, ro \f/g\ is non-zero finite. 
This is also equivalent to / ~ Ag for a non-zero constant A. 



ply that there must exist a radius below which r < 2GM,-c~ 2 . This 
would lead to a gravitational collapse to a singularity and neces- 
sarily a central black hole (or at least the relativistic effects must 
be considered for proper understanding). If the particle dark matter 
are indeed cold, the origin of the pressure support is crucial to un- 
derstand any presence of the density cusp or core. If the cold dark 
matter (CDM hereafter) halo consists of fermions, then the quan- 
tum mechanical effects, in particular, the maximum phase-space 
density is limited by the mass of the particles (i.e., df max oc mf) and 
the consequent degenerate pressure may play a role before the rel- 
ativistic effects do. In this regard, it is of special interest to study 
models with a finite central density (see Sect. 14. U which generally 
have a finite phase-space density at the centre too, if dark haloes 
are limited by the degenerate pressure of fermions. 



3.3 Dynamical properties 

Assessing the dynamical properties associated with these models 
typically requires additional assumptions e.g., regarding the ve- 
locity anisotropy etc. Some of these are straightforward: for in- 
stance, under the ergodic distribution assumption, the u nique self- 
consis tent distribution may be found by means of the Ed dingtonl 
(1916) formula although usually the tasks are more involved and 
in general the results are somewhat cumbersome. The reader is 
again referred to Appendix [B] for more systematic studies. Here 
we just note that the Hernquist ( 1990) model in this regard is the 
most special (cf. jBaes & Deionghell2002l) . for a q uite a few c om- 
plete analytically-tractable dynamical models of the Hernquist pro- 
file are known. In Figure|4]for example, we present the distribution 
function and the differential energy distribution with th e constant 
anisotropy given by = \,—\ for the dark halo of the Hernquist 
profile with a central black hole (see Appendix lB2.3l for details). 
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Figure 5. Density plots of select cored halo models, compared to the InfwI 
profile (top) and the Hernquist model (bottom). 



4 THE CUSP/CORE PROBLEM AND DENSITIES IN 
PROJECTION 

We now turn to comparisons of these models with simulated haloes, 
and consider them in the context of the core/cusp problem. 



4.1 A Family of cored profiles 

Whilst earlier N-body simulations seemed to suggest that CDM 
haloes follow the universal density profile with a central density 



cusp, emerging evidence from recent high resolution simulations 
indicate that those findings may have resulted from insufficient re- 
solving power of earlier simulations. N ewer simulations such as the 
Aquaris Project ( Springel et al. 2008) that better resolve the cen- 
tral behaviours show continu ally varyi ng logarithmic density slopes 
(cl aimed to be well fit by the Einasto profile). We have shown that 
the lEinastol profiles are an extension of the double power law to a 
member with a central core and an exponential fall-off. With a suffi- 
ciently slow transition and the consequently broad transition zone, 
continually varying logarithmic density slopes are not necessarily 
inconsistent with the presence of a central core, given the difficulty 
of simulating and resolving the central region of dark haloes. De- 
spite the lack of nominal thermal pressure support, there are some 
theoretical arguments, such as the degenerate pressure of fermion 
'gas' that argue for a core at the centre of each CDM halo. All these 
shepherd us to study the cored members amongst our model family. 

In Figure [5] we plot the density profiles of some cored double 
power-law mo dels and th e lEinastol profile compared to the iNFWl 
profile and the lHernquistl model. The models are deliberately cho- 
sen so that they closely approximate the latter models near r « r_2 
except that they are constrained to possess a core rather than a cusp 
at the centre. The finite resolution prevents any information on the 
density profile near the very centre from being available to us. For 
any practical fitting, thus there essentially exists an inner cut-off ra- 
dius below which the behaviour of the density profile is immaterial. 
Furthermore, no dark halo is in isolation in reality and so the power 
index for the asymptotic fall-off is an idealized construct. In fact, 
the fitting to the density profile should be cut off before the local 
halo density approaches the non-zero background value. However, 
there is no obvious theoretical prejudice regarding the location of 
this truncation radius relative to the transition zone of the double 
power-law profile. That is to say, the outer logarithmic slopes of 
simulated dark haloes simply correspond to that at the truncation 
radius, rather than the 'true' asymptotic slope, and there is no obvi- 
ous argument for the former to be close to the latter. With external 
factors such as tidal effects and hierarchical environments affecting 
the outer cut-off, it is thus possible for the halo to have the local log- 
arithmic slope of ss 3 at truncation as reported by many simulations 
whilst the 'true' asymptotic slope, if extrapolated, would be steeper 
(even exponential). These observations and the uncertainties in the 
fitting suggest that, given the fitting is weighted near r_ 2 , the dis- 
tinction between the conventional double power law and a cored 
profile with a sufficiently broad transition zone afforded by small p 
is nontrivial and they may even be degenerate in some degree. 

In practice, the presence of the central density cores versus cusps 
in dark haloes may be discriminated by observing y-ray signals at 
the Gal actic centre due to p air annihilations of CDM particles (but 
see also lZhao&Siik1l2005[) . In the following (Sect ED, we shall 
briefly anal yze the behavi our of the so-called astrophysical term 
J((f) (e.g.. lLapietal1l2010h in the CDM annihilation signal for the 
dark halo density profile given by our model family. We find that if 
the dark halo density profile diverges faster than r~' 12 as r — » 0, the 
astrophysical term also formally diverges as cf> — > - in particular, 
the r~ a cusp (a > ±) results in J e 0(0 _(2ar_1) ). Whilst the finite an- 
gular resolution indicates that the observed signal should be finite 
(unless the de nsity c usp is steeper than r~ 3/2 ), the signal from the 
r cusp of an NFW halo should be significantly stronger than that 
from a cored halo provided that both are normalized locally in the 
Solar neighbourhood. For example. lLapi et al] d2010h find that the 
INFWl halo produces one to two orders of magnitud e stronge r anni- 
hilation signal towards the Galactic centre than the lEinastol model 
with the resolution ranging from ss 3 deg 2 to a; 100 arcmin 2 . 
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4.2 Surface density profiles 

In terms of the comparison to the real astronomical observations, 
it is important to fold the effect of projection into our models. The 
most basic among these is the integrated mass density along the 
line of sight (i.e., the surface density profile) given by 



W) 



g(r) r dr 
Vr 2 - R 2 ' 



(7) 



This is of use for the lensing convergence k(R) of a halo or the 
surface brightness profile I(R) of a 'transparent' constant mass-to- 
light ratio stellar spheroid, both of which are proportional to ~L(R). 
If s > 1 (0 < S < 1), then equation ((7} diverges for any finite R and 
so £(/?) is not well defined. For the remaining cases that < s < 1 
(6 > 1), the result can in principle be evaluated numerically - for an 
arbitrary parameter combination, the integral in general results in 
the Fox H-function, but the particular results are in practice rather 
'useless' for the purpose of calculations. 

With the cusp/core problem in mind, the important characteristics 
of the surface density profiles of note are their limiting behaviours. 
The leading terms of Y,(R) as R — > or R — » oo may be obtained by 
analyzing equation (|7) without actually resolving the integral. The 
fundamental results in these regards are 



go) 



2 -T 

P 



(s>0) 



,,(l-»)/p 



= 0) 



hm^= B i, 

R^0Rg(R) \2 



hm 

R^°° Rg(R) 



4 



(8a) 



(8b) 



where T(u) and B(u, w) are the gamma and beta functions. Equation 
< f8ab is valid for a < 1 whereas equations l [8at are finite for a > 1 
and 6 > 1, respectively. That is to say, the surface density at R = 
is finite for < a < 1 whilst it is cusped as /J~ (a_1 ' — > oo near R = 
if a > 1 and it f alls off like Rr 16 ^ for s > (6 > 0). Please see 
AppendixlCl(also lAn & Ev"ansll2009l sect. 6) for more details. 

4.3 The Emission measure 

From the mathematical point of view, the quantity defined such that 
[g(r)] 2 rdr 



J 



f g 2 d£ = 2 f 

Jl.o.s. JR 



Vr 2 - R 2 



(9) 



shares common properties with the surface density defined in 
equation Q, corresponding to the parameter set of (2a, 26). This 
is related to the emission measure of the free-free radiation 
(Bremsstrahlung) and the so-called astrophysical term in the CDM 
annihilation signal from an external dark halo. For the dark halo of 
the Milky Way, the astrophysical term defined to be 

[g(r)] 2 rdr f r ° [g(r)] 2 rdr 

1 (O<0<f) 



p [ e (r)] 2 rdr | p 
J re (r 2 ~ 



r 

)frdr 



(r 2 - r 2 ) 



2 -J/2 



(10) 



d) 1/2 



(? < <P < n) 



with r,p = r B sin <f> is more relevant to the actual observable signal. 
Here r v is the virial radius of the Galactic dark halo, r Q is the Galac- 
tocentric distance of the Sun and <f> is the angle between the line of 
sight and the Galactic centre. 

Analogous to 11R), the general analytical expressions of J(R) or 
J((p) for the three-dimensional density profiles considered in this 
paper are expressible using the Fox H-function whilst their nu- 
merical evaluations may be achieved straightforwardly through the 



quadrature sum. The limiting behaviours of J~ as R — > are found 
similarly. In particular, for a cusped profile g e 0(r~ a ) with a > 0, 



J 6 



0(R- <2a - [) ) (a > |) 
e»(ln^') (ar=j) 
0(1) (a<|) 



(11a) 



\J(0)-J(R)\~A'R l - 2a (0<a<i), 
whilst we find JJ e 0(1) and 



\J(0)-J(R)\e 



0(R' +P ) (0<p< 1) 
0(R 2 lnR t ) (p=l) 
0(R 2 ) (p>l) 



(lib) 



for a cored profile such that \g(0) - g(r)\ ~ Ar p . As for J(if>), we 
note J((f>) ~ J7(r Q sin (p) and sin0 ~ t/> as cf> — > 0. It then follows 
that J((j)) diverges as (/r* 2 " -1 ' or ln0 - ' for a cusped halo with a > 



i or or = i. On the other hand, J((j) = 0) should be finite if the 



Galactic CDM halo is cored or cusped slower than r~ 1/2 . These last 
cases may in principle be discriminated by examining dj/d(f>, which 
formally diverges for a cusped halo as if> — > whilst it vanishes for 
a cored halo at the same limit. We leave more detailed studies of 
the projected properties of the models for the future. 
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where B(m, w) = B(l;u, w) = T(u)T(w)/T(u+w) is the beta function 
and r(u) = y(u; oo) is the gamma function. If s > | (6 < 3) on the 
other hand, then the total mass is infinite. Specifically, the enclosed 
mass diverges logarithmically if S = 3 or like r 3 ~ a if S < 3. 



APPENDIX A: ELEMENTARY INTEGRAL PROPERTIES 
Al The Enclosed mass 

The mass enclosed in a sphere of radius r is obtained by 
M, =4n | g(u)u 2 du. (Ala) 



[ I g(u) u 2 
Jo 



With our models, we should limit a < 3 since the integrand behaves 
like r 2 ~" as r — > (if otherwise, this diverges). This is directly 



related to the gravitational acceleration and the circular speed, 
vl GM, 



(Alb) 



For a < 3, p > 0, and s = 8 1 > 0, equation dAla| > results in 
AnC„ I x p 3- a 8-3 



M r 



B 



s (3-a)/p 



1 +XP 



: 6-2 r_ 2 



p 



p 

(8 - af-" 



(2 - ay-<*(6 - 2f 



i/p 



where 



a? 



Sri 



s ~~3 



2 - a I r 

\r- 2 



8-2 



(A2b) 



whilst B(z;u,w) = B,(u, w) is the incomplete beta function. The 
beta function is efficiently evaluated using series or continued frac- 
tions, which converge typically faster than integrating equation 
dA 1 ab via the quadrature sum. Some reduce to expressions involv- 
ing only elementary (algebraic and logarithmic or inverse trigono- 
metric/hyperbolic) functions, which are in general possible if any 
of the following is true: 

• (3 - a)l p or (8 - 3)1 p is a positive integer, 

• (S - a)l p is a zero or a negative integer, 

• (8 - 3)/p is an integer and (3 - a)jp is a positive rational, or 

• both 2(3 - a)lp and 2(6 — 3)1 p are integers. 

Equation dAlat for s = on the other hand results in 



■y 



3-a 



(A3a) 



l-A/,0 - Qs'sP 



where a < 3, p > and 
* 2 - a I r x '' 

r_2 



y 



(A3b) 



pn p 

Here, y{u;z) is the lower incomplete gamma function. This 'sim- 
plifies' to elementary (up to exponential) functions if (3 - a)/ p is a 
positive integer. It is also recast to formulae using the error function 
and elementary functions for half-integer values of (3 - a) I p. 

The behaviour of the integrand, r 2 ^ 6 as r — > oo of equation dAlat 
indicates that the finite total mass is defined for 6 > 3 (including 
the s = case). In particular, 



4nC M J3_ 



P 

4nC M ,o 



a 8-3 



P 

3-a 



P 



(0<s< |) 
(* = 0) 



(A4) 



A2 The Gravitational potential 

The potential of an isolated system up to an additive constant is 
found by integrating Vi/f = -g. Under spherical symmetry, 



# = GMr = _i_/^M + 47tG r 
dr r 2 dr\rl 



(A5) 



where we have used dM r /dr = 4ngr 2 . Then the potential behaves 
asymptotically like r -1 , lnr, r~ <6 ~ 2> and In r for 8 > 3, 8 = 3, 
8 < 3 but 8 + 2, and 8 = 2, respectively. If 8 > 2, the potential at 
infinity is not divergent and any particle with enough kinetic energy 
can escape to an unbound state. By contrast, if <5 < 2, then if/ — > oo 
as r —> oo, and thus every tracer particle in the model with i' > | is 
bound to the system and escape is formally impossible. 
If < s < j (<5 > 2), it is customary to set (/'(oo) = and thus 



(A2a) ^( r ) - _ AnG J u du 



For s = S-' > 0, 
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AnG 
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& 3 Si 



1 + x p p p 



+ B 



C M 2 s, lh , 

8^>r s &Ss 

whilst, for s = 0, 



AnG 



C, 



ip,0 
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Cm. a 
p l/p r s 



1 



«3 



; Q-2 r_ 



y" +r 



(8 - a) 6 



1 82 &2 
1 + XP ' p ' p , 

Up 



(A6) 



(A7) 



a 2 



(A8) 



■e s r 2 P CnJ " 



■ Q-2 r_ 



& 2 /p 



S 2 /p 



where F(u;z) is the upper incomplete gamma function. We also 
use short-hand notation, &„ = (n- a) and 8„ = (8 - n) for an inte- 
ger n. These again reduce to expressions involving only elementary 
functions if M r is expressible in such a way and the parameter com- 
binations 82 1 p and ch/p satisfy criteria similar to those applied for 
83/ p and & 3 /p. Some specific examples are found in Appendix ID1 
The integral in equation iA6i for 8 < 2 is divergent, and so equa- 
tion dA7t for such are meaningless (although they are technically 
well-defined unless <5 2 / p is a negative integer). The potential for 
8 < 2 however can still be defined with an alternative zeropoint. 
Since M r e &(r 3 ~ a ) as r — > 0, we have if/ e &(r 2 ~ a ) if a t 2 or 
<9(lnr~') for a = 2 in the same limit. Hence, provided that a < 2 
(note that we have generally restricted to be a < 8), the alternative 
reference for the potential <^(0) = is well-defined and 



, GM,. 
<M>)-^<= — +4ttG 

For s > 0, this results in 



I g(u)u 
Jo 



du. 



AnG 



Cfi 
P 



1 



B 



x" _ q 3 Si 
1 + xP ' p ' p , 



(A9) 



(A10) 



, 1 + x'' ' p ' p ( 

If a < 2 < 8, both this and equation dA7t are well-defined, 
but they differ by a constant (independent of r) - note B(\ — 
z\u,w) + B(z;w,u) = B(u,w) for u, w > 0. This difference cor- 
responds to the depth of the central potential well, \if/ m - ifi \. This 
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and -y/2 \ijfZ - >Po\ at" e also the maximum specific energy and speed 
of any bound particle in the system. 

For s = 0, equation dA8t is always (provided that a < 3) well- 
defined. That is to say, since the total mass is finite, it is always 
possible to set i//(oo) = 0. The central potential well depth is then 
the limit of equation {A8) as r — > 0, which is finite if a < 2. Using 
this and y(u; z) + f(u; z) = r(«) for u > 0, it is possible to derive an 
alternative expression for equation d A8t for a < 2, that is 

- if _ Q.o 
P 



4nG 



1 &3 , 

-y — ,y 

y \ p 



a 2 , 

■y — 

p 



(All) 



of which the zero point reference is essentially at the centre. 



This inclu des the y-sphere for which p = 1 and 6 = 4, the model 
studied bv lEvans & Ar] fe005t) for which (a, 6) = (2 - p, p + 3), 
and the analytical solutions of iDehnen & M cLaughlin! d2005T) for 
which (a, 6) = (1 + A> - f.P + 3). Note B(z;u, 1) = z"/u, and so 
<x 2 for this last case results in a double power law as the density 
profile. The resulting model with p = 2(2 - 8 )(e - 2)/(s + 6) 
exhibits radial powe r-law behaviour, p/crf oc r~ n w here r\ = as/2 + 
a-Bge as noted 2 ! bv lDehnen & McLaughli Equation l lB6t 

also indicates that cr 2 ~ cr 2 /(3 - 2BJ) ~ v 2 /(p + 4- 2Boo) as r — > oo 
whereas, as r — » 0, we find that a], e 0(r°' _2 ' ?o ), &(r 2 ~ a \nr~ [ ), and 
0(r 2_a ) for 1 +y6o > cf, 1 +/3o = a, and 1 +0p < a respectively. T hese 
are consistent with behaviours deduced by I An & Evans! ( 2006). 



APPENDIX B: TRACTABLE DYNAMICAL PROPERTIES 
Bl Velocity dispersions 

The velocity dispersion profile of a spherical system is found by 
integrating the spherical steady-state Jeans equation, 



1 d( ec rf.) 



GM r 



g dr r r 2 

Equation ( IB It may be integrated to be 



0. 



<t 2 M) 



du 



2f3(u) 



(Bl) 



(B2) 



given the boundary condition lim,^ m gcr 2 = 0, once g(r) and M r , 
as well as the so-called velocity anisotropy parameter 



m = i - 



2a 2 



(B3) 



are specified. For a monotonicall y-varying anisotropy par ameter, 
the parametrization introduced bv lBaes & Van He"se] J2007t) . 

Poor 2 " + p r: q d log[r 2 »(r 2 « + ,*)<&»-/%)/»] 



r 2 i + rl" dlogr 2 (B4) 

may be useful to represent its behaviour. Then, equation dB2t re- 
sults in a simple integral quadrature for cr 2 {r), namely, 



h(r). 



i r°° gm u 

— — T-h(u)du; 

Mlq g{r). 



K 



(B5a) 



(B5b) 



which can be evaluated numerically if M r is tabulated. Since h 6 
0{r 2l3a '^ 6 ) as r —> oo, th is is convergent pro vided that 1 +/3oo < S (cf., 
Hansen 20o3, see also I An & Evanslfeood) . Given /? < 1, equation 
JB5ab is therefore we l l defin ed for 5 > 2. 
lAn & Evansl |2006, 2009) have analyzed the leading term be- 
haviours of a 2 in the limit of r — > and r — > oo. The complete 
analytic expressions for arbitrary parameter combinations are yet 
possible only with higher transcendental functions, which are more 
or less impractical. For some parameter combinations however, re- 
ductions to simpler functions are a vailable. If (3 is a linear func- 
tion of the logarithmic density slope dHansen & M oore 2006), then 
r a = a and 2q = p, and so h(r) in equation dB5bt also takes the form 
of a double power law for s > 0. If M, is additionally reducible to 
the sum of such factors (e.g., 63/p or a 3 /p is a positive integer), 
then a 2 is expressible using the sum of incomplete beta functions. 
The simplest example is the 6 = p + 3 case, for which 



4nGC, 
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x"- 2llo (l+x") 



y 2(fe,-l0)-S 3 



x B 



1 2(2 -/U , 2(l+/? -or) 



1 +xp 



+ 1, 
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(B6) 



B2 Constant-anisotropy distribution functions 

Consider the phase-space distribution given by 



r(B, l 2 ) - 



7r 3 ' 2 r(l -^L 2 " 
6(^ 2 )g(£) 

V27T 3 / 2 



08 < 1) 



03= 1) 



(B7) 



where fi and L are the specific binding energy and the magnitude 
of the specific angular momentum and 6(w) is the Dirac delta. This 
builds a spherical system with the anisotropy parameter at all radii 
given by the parameter value B. The function g(£) is related to the 
tracer energy distribution via 



d£ _ i^-si^epF-s) 

d£ " r(§-j8)rV 
dM _ N(8)g(S) 

dS ~ ra-0) 



N(£) = 4n 



(B8) 



\V - r^dr, 



where @(w) is the Heaviside unit step function and r E is the farthest 
distance to which a particle with the binding energy £ can move, 
that is, v P(rg) = £. Here, dg/d£ and dM/d£ are the local and global 
differential energy distributions and ./V(£) is the marginalized den- 
sity of states. The function *?(?•) defined to be 



^out) - f (r) 



if r ou , is finite 

if ''out = 00 and I'M 00 )! < 00 
if '"out = 00 and (/'(oo) — > 00 



(B9a) 



is the positive potential function where r out is the boundary radius. 
The specific binding energy is given by £ = - | u 2 where t) is the 
speed of the tracer with its lower bound for bound particles being 







for an escapable system 



-00 for an infinitely-extended inescapable system 



(B9b) 



The undetermined function g(£) is uniquely specified by the den- 
sity profile. That is, integrating equation dB7t over the accessible 
velocity space at a fixed location results in the augmented density 
hiW) = r^g as a function of ¥, The function g(£) then may be 



2 The complete specification of the model involves 7 parameters, but the 
choice of is independent from the others and S = p + 3 indicates 
that specifying either p or 6 also fixes the other. This effectively leaves 5 
parameters and 3 algebraic constraints, and thus specifying two amongst 
(a, p,j8o, € , T]) completely determines this model. 
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determined by its inversion, the lEddingtonMCuddeford formula, 



g(S) = 



1 



r(l - e) dfi 
1 



= T(l - e) 

h ( "\S) 



if 

J s ~W- 



(S - >F) e 
WdY 



(fi - £) £ 



(0 < 6 < 1) 
(6 = 0) 



(BIO) 



where p = |_§ - /3J and e = (| -/}) — ft are the integer floor and 
the fractional part of | - /?. Equation l IBlOt indicates that, given 
the inverse function r(if/) = of if/ = if/(r) and the augmented 

density h(?¥) = [r(- > P)] 2/? £)[r(- , P)], the distribution function (df) of 
the form of equation ( |B7t is available as an integral quadrature. If 
/3 is a half-integ er (i.e., B = ^, — 4, . . . ), the formula results in pure 
differentiation iCuddefor j[l99ll ; jEvans & Anll2006h . 



B2. 1 the (a, S) = (2- p, 3 + p) case 

The simplest example is the (a, S) = (2 - p, 3 + p) case, for which 
ga 2fi <t> p+3 - 2fi _ ¥ (p + 1)T 



(1 _ (D;>)2(i-«/;.-l 



l^oo - tAol 4ixGga 2 



(Blla) 



Here, \tf/„ - (fr l = GM tot /a = 4nGga 2 /(l + p) is the depth of the 
central potential (note < <£> < 1). The df is then found with 



*(£) 



gcr 
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-</fo\ 



r(l-e)J (®-0 f de' 



)i, + 1 



gp+3-20 



(1 _ Qpyxi-pj/p-i 



(Bllb) 



If 2(1 — B) > p > 0, the df is non-negative for all accessible bind - 
mg energy interval, < S < l^-^ol (cf.. lAn & Evan3l2006h . 
Equation JBllbb with arbitrary p > and /? < 1 results in the Fox 
H-function, which is equivalent a power series of O, namely, 



(8 p ) k r(pk + p + 4-28) 
kW{pk + p + \-p) 



(B12) 



where /3 P = 2(1 — B)lp - 1 > and (a) k is the Pochhammer symbol. 
This re duces to a h ypergeometric functi on if p is rational. The resu lt 
for the Hernquist model was derived bv lBaes & Deionghel (2002), 



2 Fi(l -2/3,5-2,8; 2 -#*), 



r(5 - 2/3) 5 

r(2-^) 

whilst for the I S chusteJIPlummerl sphere, this becomes 



(B13) 



§(<*>) = *KF 2 (-ft3- A | -A ^, ^;0 2 ). (B14) 

If /? is a half-integer, both are expressible using only algebraic func- 
tions whereas equation dB13b with an integer /3 reduces to an ex- 
pression inv olving up t o inverse trigonometric functions, e.g., for 
the isotropic lHernquistl model, 



2^n 



3 arcsinV* V* ( 1 - 20)(3 + 80 - 8<D 2 ) 



(!-<»)§ 



(1 - ®) 2 



(B15) 



If 2(1 - /3) = p, equation JBllbt indicates that g(O) oc <bP + i-P = 
^)<3p+i)/2 Equation ( IB 121 is still valid because the series then ter- 
minates for k > 1. The correspond ing df, T(S,L 2 ) oc L p - 2 6 (3p+1)/2 
is the hypervirial model studied bv lEvans & AnU2005l) . 



B2.2 the distribution function as a function ofrs 

Even if an analytic expression for the inverse function of the poten- 
tial is not available, equation dB lOt can be used to write down g(£) 
implicitly using r$ = ¥ (S) = if/~ '(-£). This is achieved through 
change of variables for the integral in equation dB 10b . 

1 d^ +1 f 00 _ GM r g(r)dr 
" ' m-e)d&^J rp r 2 

(B16a) 
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ni-e)Mr E dr e 
r 2 d 



aa-»[fi - ¥( r )Y 
s M r ^{r)dr 



[r^gir)}. 



(B16b) 



^ \GM,-dr) 
Here *¥ = -iff and dtfj/dr = GM,jr. In addition, v P(r fi ) = £ and so 



dfi 
dr E 



dip 
"d7 



GMy 



d 

dfi 



GM r Are, 



The integral in equation dB16at for a positive integer n = | - 
collapses thanks to the fundamental theorem of calculus and so 

g(£) = (-ly^foj). (B17a) 

For/? = i (n = 1) and/3 = — A (n = 2), this reduces to 

dlogg^ 



rg 



(GM r ) 2 



dlogr 
d 2 log g I Anr^g 



dlogr 



dlogg 

M r d log r 



1 - 



dlogp^ 
d log r i 



The df and the differential energy distribution are given by 



and 



dM 
dfi 



2nr\g(&) for/? : 



(B17b) 



(B17c) 
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Jo 



GM r rdr 



T(S,L 2 )=—g(S); 
dM n- , 

— = 47tg(£) J o CP - fi)r 3 dr = Ttg(fi) 

for /3 = - j . Writing down the correspondin g df for our models 
using equations dB 1 7b is trivial if tedious. See |E vans & Arj d2006h 
for simple cases such as (p, 6) = (1,4) and (p, a) = (1, 1). 



B2.3 a central black hole 

The outlined procedure is valid even if there exists a central point 
mass a la a black hole. This does not contribute the local density 
except at the centre, and so it is simply incorporated by M r = 
^r.haio + M, and tf/ = i/vio ~ GM.r~ [ , that is, the mass of the central 
black hole M. is an integration const ant. With n on-zero M. > 0, 
the potential if/(r) is invertible for the iHernquisll model and those 
correspond ing to (a, 6, p) = (|,4, 1) and (1, f , !)• If/ 3 = \, the 



Hernquist model with a central black hole (Ciot til 19961) is built by 
r(£L 2 )= 3 «; 2 (l-«;) 2 dM 3a w 4 



4n 3 aGL w 2 + . 



dfi G w 2 + .4 



\-g-jZ+ ^(\-£-Jf) 2 + 



(B18a) 



For (a, 6,p,B) = (|,4,l,i), we find teaes & Deionghell2004h 

3 (1 - io 2 ) 2 (l + w 2 ) dM 3a iu 3 (l + w 2 ) 



T(S,L 2 ). 



32n 3 aGL w(w + yi 



dfi 8G w + . 



i + V(i + + sj( 



2 + g + . 



(B18b) 
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whilst the model with (a, S,p,fi) = (1, |, 1, |) is given by 

r(£ T2, = z(l-z) 2 . 

(27t) 3 aGL4(l- Vz) 2 + Vz^' 
dM _ 3a (1-z) 4 
d£ ~ 2C Z [4(1 - Vz) 2 + Vz.^]' 



(B18c) 



_ 8(4 - <r> - (16 - <r - + 8 ^{4-~¥fTJ^M 

W ~ (8 - S - ^T) 2 

with z = 1 - t«. For all cases, S and are the binding energy and 
black hole mass normalized such that 
(3 - a)6 M (3 - a)M. 



S 



AnGgd 2 ' 



4nga 3 



Calculations for p = —j are sti ll trivial a lbeit messier. We only 
provide the explicit result for the Hernquist model, 

2 _ L(l-ai) 3 u> 2 (3 + 4w + 3W 2 ) + J£(\ + 3w + 6w 2 ) 
,L ] ~ 2n 3 a 2 G 2 M tot (w 2 + J(f 

dM _ a w 2 (3 + 4w + 3w 2 ) + .y£(l + 3w + 6w 2 ) 
dfi " 3G (w 2 + J(f 

x \u)(Yl - 30w + 22iu 2 - 3u?) + .£u? + 12(1 - w) 3 ln(l - w)] , 



where w is as in equation dB 18at . and $ = a£/(GM tol ) and = 
M./M m where M tot is the total halo mass (sans the black hole). 



APPENDIX C: LIMITING BEHAVIOURS OF X(R) 

For an arbitrary parameter combination, equation {7} results in 

(1-6,1) 



g s r s \Rl 



x 



V^tt2,1 

r(b) " 2 2 



V^H 2 ;» 



la £i 
l 2 . 2 ) 



{0, !},{"=*, 



(!.§] 



{0,l},{Hji,f} 



(CI) 



1/R\ P 

where b = (6 - a)/ p. The first case is for < s = 5~ A < 1 whilst 
the second case corresponds to s - 0. For a rational value of p, this 
Fox H-function reduces to the Meijer G-function, and subsequently 
it is possible to write down using the hypergeometric functions, if 
desired. If p = 2 or p = 1, this reduces to simpler special functions 
or elementary functions for some parameter combinations. More 
detailed studies of these types will be explored elsewhere. 
We also find that if 1 < a < 6, then 

ICR) 



e*r s \RI 



V^r(H=i) 



r(f) 



+ 0(g) 



(R->0) 



(C2) 



where the reminder term g(R) behaves like 

Rf (p<a-l) 
g(R) = • R"lnR [ (p = a - 1) . 

R"- 1 (p>a-l) 
If a = 1 < 6, then I(i?) diverges logarithmically like 
I(if) , ,2r s ^ y-F 



2f? s r s 



+ 0(g), (tf-»0) 



(C3) 



\lnS- 
Unp 



F( S -f) (»>0) 
(ff = 0) 



where F(z) = dlnT(z)/dz is the digamma function and y 
0.5772 ... is the Euler-Mascheroni constant with the reminder 



g(R) 



RP {0<p<2 K p±V) 

R 2 lnR i (p=\ V p = 2) 
R 2 (p > 2) 



For 0<o-<l<<5on the other hand, it can also be shown that 
X(R)-I(0) 2V7tr(2fi) ^i-" 



Q s r s 



(1 -a)T(f)\rJ 



(R -» 0) (C4) 



A 2 In A 1 
tf 2 



(p < a + 1) 
(p = a + 1) . 
(p > a- + 1) 



In other words, l.(R) for < a- < 1 decreases outwardly from 
the central line of sight like |I(0) - I(R)| ~ AR l ~ a and possesses 
a 'spiky core', i.e., dl/dR diverges as R — > 0. With !(/?) corre- 
sponding to a cored profile with a = 0, we find that it exhibits a 
'flat-topped core' at the centre as 



1(0) - 1(g) 

s r s 



ln(^)-H + i](^) 2 +0(^ 4 ln^ 1 ) 

Ui + P)r(i-f) U) 



G : 



<5 1/p r(^) 



r(^) 



(.5 > 0) 

(s = 0) 



H = 



Y 



In 5 (,s>0) 
(j = 0) 



0»>i) 
(p = i) , 

(P<i) 

(C5) 



where //„ = F(m + 1) + y is the harmonic number. That is, given 
the cored 3-d density profile \g(0) - g(r)\ ~ Ar p as r — > 0, the cor- 
responding surface density as R — > behaves as |X(0) - I(R)| ~ 
A / / jmin(i +P .2) unless p = i > for w hich it does as |I(0) - Z(R)\ ~ 
A'R 2 In R l . The order of the reminder term is given by 



g{R) 



The limiting behaviour for R — > oo on the other hand is easier to 
derive because the Fox H-function in equation JClb for b > and 
s t reduces to a power series, 



R 4 


(p > 3 V p = 2) 


R 4 \nR- 1 


(p = 3 v p = 1) 




(1 < p < 3 A p + 2) 


R 2 


(\<P< 1) 


R 2 lnR- [ 


(P= |) 


R 2 P+ i 


(0<p< i) 



I(R) = (r.\ a ~ l Vn y r(-^- 
p s r s \W/ X 5 - 0, 4^ rY£+ffi 



)(-D"(t £ )„ 



r(2^) 



(C6) 



which converges for |X| > 1. Here X = S~ [/p (R/r s ). This is consis- 
tent with equation d8bt . For the s = case however, equation J8M 
indicates that 1(A) falls off faster than Rg(R) and we instead find 
the following asymptotic expansion (R — > oo) 



Z(R) 



/>/2 



1 - 



4a + 3p - 6 /r s 



(C7) 



yl2jxRg(R) Vs 

That is to say, we find that S(A) ~ yl2jixR l -'i g(R) for a 3-d density 
profile behaving as g(r) ~ Ar~"e~ Br '' . 
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density profile fitting functions 1 1 



Table 1. Potential-density pairs for a = 2 - mp and 5 = 3 + np. 



m n Cq- 1 -(47tGC) - V 



1 1 , 2 -P(\ +xi') 1+T ' 
1 2 



2 



1 +p (1 +jcP)'/p 

r2-»(l+^) 3+ 5 £ !_fp + _!_) 

(1 + p)(l +2p) (1 + xP)Vp V l+xP) 

i ^(i + ^) 3 4 ^ L_L + ^L) 

(1 + p)(l +2/)) (1 +xPflP V l+xP) 
1 3 + - ■ - g — — — ' , /n [2p 2 + - 



(1 + p)(l + 2p)(l + 3p) (1 + ^p)Vp L * 1+jcP 

2 2 r^n+x") 4 ^ ^ ! \p+ X " 

3 i ^ ( i +x pfH a - l_ r 2p2+ i^ + ( i+ ^; 



(1 + + 2p)(l + 3p) (1 + xP) 1 /? L 1 + (1 + ,v' J ) 2 



Table 2. Potential-density pairs for or = 3 - mp and 5 = 2 + np. 



m n Cq 



-(4ttGC)-> 



1 1 r 3 -''(l +.y'') 2 ~? 

1 2 r 3 -''(l + x p ) 3 ~r 

2 1 r 3 - 2 P(l +xP) 3 ~? 
2 2 + x p f~? 



(1 + .v p ) 1/ '' /l 1\ 
r \ r a / 



(l-p)(l-2p) 

n 2 P 



I 1 \(1 +xP) [/ P p-l p 

\l+xP P ) r + r + a 
/ xP \ (1 +x p f lp p p-l 
\ 1 + xP I r r a 

U xP 1 (1 +xP) 1, p p p 

(l-2p)(l-3p)\i(l+jfP) 2 P \ r + r + a 



(l-p)(l-2p) 
,P-p 



If the potential results in § overall, it reduces to an expression involving the logarithmic 
function by the continuous limit. See Tabs.|6](for p = 1) and|7](for p = j). 



APPENDIX D: ANALYTIC POTENTIAL-DENSITY PAIRS 

Thanks to the properties of the beta and gamma functions, any one 
of the followings in the list is the sufficient condition for the poten- 
tial for our model family to be written down analytically, namely, 

(i) (a, 5) = (2 - mp, 3 + np) or (a, 8) = (3 - mp, 2 + np) where 
m and n are positive integers (Tabs.[T]and[2]l, 

(ii) p~' is a positive integer, and cb/p or 63/ p is a positive inte- 
ger (Tabs.[3]and[4j, or 

(iii) all of 21 p, 2a/ p, and 2r5/p are integers (Tabs. l5l7l l. 

Tables [T] and [2] gi ve examples of the case (i), which include the 
models studied by lEvans~& Aril (2005) corresponding to (a, S) = 
(2 - p, 3 + p). The examples for the case (ii) a re found in Ta bles [3] 
and[j] These include t he so-called y-sphere of lDehner] d 19931) and 
Tremaine et al. ( 1994) for which 63 = p = 1 and the generalized 
NFW profile studied bv lEvans & Anl (2006) for which & 2 = p = I. 
This last family is extended to include the model with i' = 0, i.e., the 
potential-density pair of if/ oc (e~'' r — l)/r and g oc r~ l e~ w . Finally, 
in Tables 15171 we list some examples of the case (iii) with p = 2 
(Tab. [5}, p = 1 (Tab. |6j, and p = | (Tab. |7J, respectively. Some 
examples are in fact redundant as they also belong to one or both 
of the first two cases. In Table [7] we only list those with both 2a 
and 26 are integers, but if both 4a and 46 are integers with p = \, 
the resulting potential is still analytic. However, if 2a or 26 is a 
half-integer in this case, the resulting expressions involve inverse 
trigonometric/hyperbolic functi ons. 

The potential generated by the lEinastal profile with an integer in- 
dex n = p~ l is another example of the case (ii) that is expressible 
using only elementary functions. For the s = cases, that p = n~ [ 
and a = 2 - m/n where both m and n are positive integers is a 



sufficient condition for the existence of an expression for the po- 
tential written using only elementary functions, which corresponds 
to the case (ii). On the other hand, if both 2/p and 2a 2 /p are pos- 
itive integers and s = 0, the potential can be written down as an 
expression involving the error integral. Some examples are found 
in Tables[5l7l 

Examination of these tables reveals a symmetry between the po- 
tentials for the models with (a, 6) «-> (5 - 6, 5 - a). This may be ex- 
plained as follows: Let us suppose that the potential function if/(r) is 
generated by the double power-law profile so that it is the solution 
of the spherical Poisson equation, 



id, 
r 2 dr 1 



dr 



4nGga" 



+ x p) (s - a)l P 



(Dl) 



where ga" = g s r" and x = r/a. On the other hand, the same Lapla- 
cian operator acting on the function <p(r) = r~ [ if/(a 2 /r) yields 



1 d 

7- d7 



1 d 

w 2 dw 



1 



» dif/(w) 
dw 
4nGga 5 - s 



w=a- 1 r 



(D2) 



a r 5 ~ 6 (l + x P)^~ a ^ p ' 

In other words, if if/(r) is the potential corresponding to the dou- 
ble power-law profile with the parameter set of (a, 6, p), then 
the potential for the parameter set of (5 - 6, 5 - a, p) is given 
by (r/ay l tf/(a 2 /r). This also indicates that the potential for the 
model with a + 8 = 5 (which includ es the hyp ervirial models of 
Evans & Aij l 2005l encompassing the Hernquist and the ISchusteil - 



Plummej sphere) satisfies the relation ai{/(a 2 /r) = ri//(r) for all r 



with the properly chosen scale length a. 
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Table 3. Potential-density pairs for p = l/n and 5 = 3 + m/n. 



S p Cq 



-(47tGC) _ V 



4 1 r'^l+.T) 4 -" 

5 1 r ff (l + xf- a 



1 r ff (l + V*) 2 * 4- ^ 



(2 - or)(3 - a) 
„2-a 



(l+j 



(2 - a)(3 - a)(4 - or) 
4a 2 -" 



2* + 4 - a 



(4 - 2a)(5 - 2o0(6 - 2a) 



(1 + x) 3 - ff 



2 I 

2 2 



r ff (l + V*) 9 



(-"(1 + Vjt) 2(5 ~ ff) 



(4 - 2or)(5 - 2a)(6 - 2a)(7 - 2a) 

Ha 2 -" 



(1 + Vjc) 5 " 2ff 

_ 3x + 6(3 - a)V* + (5 - 2a)(7 - 2a) 



(i + v^) 2(3 - ff) 



(4 - 2a) 5 



, 6jfJ + 6(7 - 2a)jc + (4 - or)(29 - 10a)VJ + (4 - a)(5 - 2a)(7 - 2or) 2 _ 

6 = = ; X 

(i + V*) 7-20 " 



2-n 



(4 - 2a) 6 



30- 



30a- 2 + 60(4 - a)aS + 9(9 - 2ar)(10 - 3a)x + 2(4 - a)(9 - 2a)(20 - la)ypc + 8(4 - <*)(§ - a) 3 2 _ a 

(i + vi) 2 < 4 -«) * 



Here, (a)„ = (a)* is the Pochhammer symbol. The potentials corresponding to (a,p) = (2, 1), (2, j), and (|, i) are found by the continuous limits, 
which result in the expressions involving the logarithmic function. The specific results are provided in Tabs.[6](p = 1) and[7](p = 5). 



Table 4. Potential-density pairs for p = l/n and a = 2 — m/n. 



a p Cq 



-(47tGC)"V 



3 1 
2 2 



1 1 r(l+xf- 1 
1 (l+xf 

A(i + v*) 2S ~ 3 

r(l + VA : ) 2( ' 5 - 1) 
Vr(l + V*) 2 " 5-1 

(1 + v^) M 



1 j_ 

2 2 







(5-3X5-2) r 



1 



(1 + x) 6 - 3 



2- 



2 + (<5 - 1)a 



(5-3)(5-2)(5- 1) r|~ (1+a)" 5 - 2 

4a 3/2 1 [" 1 + (25 - 5)V* 

(25-6)(25-5)(2<S-4)7[ 1 " (1 + VD 2 " 5 " 5 

4a 2 1 

(25 - 6X25 - 5X25 - 4)(25 - 3) r[" ( 1 + V*) 2 '*" 2 ' 



3 + 6(5 - 2)V* + (25 - 3)(25 - 5)x 



Sa 5 ' 2 1 
(25 - 6) 5 r 

8a 3 1 
(25 -6) 6 7 



6 + 6(25 - 3)y[x + (5 - 1)(105 - 2l)x + (5 - 1)(25 - 3)(25 - 5)jd 



30- 



(1 + v^) 2M 

30 + 60(5- 1)V* + 9(25- 1)(35 - 5)x + 2(5 - 1)(25- 1)(75- 15)jc§ + 8(<$~ 1)(5- |) 3 a 



(1 +Va : ) 2( ' 5 -' ) 



The potentials corresponding to the zero denominator for ifi result in a logarithmic expression through the limiting process if 5 > 2. The results 
are given in Tabs.|6](p = 1) and[7](p = j). The potential at the limit of 5 — > 00 (j = 0) are also found in the same tables. 

Table 5. Analytic potentials with p = 2. Here, we list -(4nGCa 2 ~ a Y V f° r 
5 < 00 whilst the corresponding density is given by Cg- 1 = r a (\ + x 2 ) (s - a)/2 . 



p=2 a=0 



a = 1 



a = 2 



arsinh x ln(l + x 2 ) , 
5 = 3 Y arctan - 



arsinh x 



2x 



5 = 4 
5 = 5 



1 +x- VT 



2x 
1 



+ arsinh - 



■ + ita(l + i) 



3VT 



arsinh ■ 



Vtt ,. r 1 - e 

s = erf y erfc y H 

4y 2 2y 



^erfy+ |£iCv 2 ) 

2y 



For s = 0, we list —{AnGCa 2 ~ a )~ l ij/ whose corresponding density is given by 
C(f l = r" exp(y 2 ). These actually do not result in potentials that are express- 
ible by elementary functions, but the expressions involve the error integral or 
the exponential E[ -integral (for a = 2). 
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Table 6. Analytic potentials with p = 1, The normalization is the same as Tab. [3] 



p = 1 a = 



ff = 1 



a = 2 



5 = 3 



6 = 4 



6 = 5 



16 



4 + 3a \ 1 



V4V1 + x 
ln(l+x) 1 



2(1 + x) 



6 = 2 is 



4 + 5x 



4(1 +x) 3 ' 2 



1 + 2.V 



6(1 + .v) 2 



u — 7 ins 



4 + 7x 



4(1 + x) 5 ' 2 
1 + 3x + x 2 



12(1 + x) 



■13 



1 arctan VI \ 1 
3 — — + arctan ■ 



VI 



' v*(* + 1) 



V3E 

2 arsinh V* 



/ 1 3 arctan -Jx x + 3 
r arctan — — H 



T5 ' 



+ i)V? 



/ 1 arctan -Jx\ 3x 2 + 2x + 3 

arctan — — H — - 

\ V* * / 24VI(x+l) 2 

(4x + 7) x 3/2 



4(x + l) 5 ' 2 



4 V-i+ 1-1 



ln(l + x) 



I 1 



1 



1 1 1 



ll+xl x 



2(1 +x) 



1 



(1 + x) 3 ' 2 



2 + x 
6(1 + a-) 2 



1 arctan V* 1 
2 — — h arctan 



, , / 1 + x arsinh VI 



arctan "Jx 1 1 
v arctan — — 

x VI V* 



, / arctan VI 1 
i h 3 arctan — — 



1 + 3a 



VI (1 + x)ypc 



(Ax + 5)VI 



4(x + l) 3 ' 2 



(ix+T-i i 

2 h arsinh — 

I x VI 

/ 1 \ ln(l + x) 
ln(l + -] + — 



x' 



X 



l-yl+x l 
2 1 h arsinh — 

x vi 



In 1 + ■ 



, i 1 + 3x \ 1 l 
I 1 ■ - + 2 arsinh — 

VT+1 j VI 



In 1 + i) - 



1 



2(1 + x) 



4 

VI 

2 i^±I_ 1 + ^W!j 

1 arctan VI 1 \ 

2 1 — — H ai'ctan — 

VI * VI/ 



1 l \ arctan -yl 

3 1 — — - ai'ctan — + 

VI VI/ x 



16 



4x + 3 



4Va(a+ 1) 



;[l-d + |)ei ^[l-d^)erfVy]-^ 



1 -e - 



VS[l-(l-^)erfVy]-^ 



Table 7. Analytic potentials with p = A. The normalization is again the same as Tab.|3]and|6] 



5 



-0 



j a = 1 a=| a- = 2 a = | 



60+ 153VI+ 119x + 24xi 10T 2(12 + 19V* + 6a) 8^+ 2(3 + 2yfx) 6t+ 4 4( + 2 2£+ 4 



6VI(1 + VI) 3 * 3VI(1 + VI) 2 * VI(1 + VI) * yfx x y[x x VI 

2f+ 60 + 207 VI + 248a + 107x2 2f+ 12 + 29VI+19x 2l + 2(3 + 4i/I) 2f+ 2 2£ + 2 2 | 21+ 

a 30VI(1 + VI) 4 * 6VI(1+VI) 3 x 3VI(l+VI) 2 * VI(1+VI) a vi VI * 



«=! 




c = 2 




4 

VI" 


X 


2 

— +2T 

VI 


2£+ 

A 


2t + 


2 


2£ + 

2r + — 

A 


2 


X 


VI(i + VD 


"VI 



1 +5VI+ 10a + 10x1 1+4VI + 6a 1+3VS 2 2 4 

v v v zr — - 4r 



7 

30(1+ VI) 5 15(1+ VI) 4 6(1 + VI) 3 3(1 + VI) 2 1 + V* VI 

3 + 18VI + 45A + 60Ai + 10x 2 1 + 5 VI + 10x + 2xJ 1+4VI + X 3 + VI _ 2(4 + 3 VI) 2(2 + 3 VI) 

2t 



210(1 + VI) 6 30(1 + VI) 5 10(1 + a/I) 4 6(1+a/I) 3 ~ 3(1+a/I) 2 a/I(1 + VI) 

3 + 2lVI + 63x+ 105a-i + 35x 2 +5xi 2(1 + 6y/x + 15a + 6xi + x 2 ) 2 + 10VI + 5x + xJ 6 + 4-^/1 + x _ 19 + 29VI+12x 2(6 + 19V* + 12x) 

2t — — — — — — — — ot 



420(1 + VI) 7 105(1+ VI) 6 30(1 + VI) 5 15(1 + VI) 4 6(1 + VI) 3 3VI(1+VI) 2 

1 + 8VI + 28x + 56xi + 28x 2 + 8xi + x 3 5 + 35VI + 105x + 63x? + 21x 2 + 3a? 10 + 60VI + 45a + 18xJ + 3a 2 10 + lOVI + 5a + a? 107 + 248VI + 207x + 60x2 24 + 119VI+ 153a - 

252(1 + VI) 8 420(1 + VI) 7 210(1+ VI) 6 30(1 + VI) 5 " 30(1 + VI) 4 6VI(1+VI) 3 



+ y[l-(l + ^S^H f[l-(l + ^ + f)e^] l[l-(l+VI)e^] 



Here, t = ln(l + VI) and C = ln(l + -J=). 



